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Abstract 

The phenotypic equilibrium, i.e. heterogeneous population of cancer cells tend¬ 
ing to a fixed equilibrium of phenotypic proportions, has received much attention in 
cancer biology very recently. In previous literature, some theoretical models were 
used to predict the experimental phenomena of the phenotypic equilibrium, which 
were often explained by different concepts of stabilities of the models. Here we 
present a stochastic multi-phenotype branching model by integrating conventional 
cellular hierarchy with phenotypic plasticity mechanisms of cancer cells. Based on 
our model, it is shown that: (i) our model can serve as a framework to unify the 
previous models for the phenotypic equilibrium, and then harmonizes the different 
kinds of average-level stabilities proposed in these models; and (ii) path-wise conver¬ 
gence of our model provides a deeper understanding to the phenotypic equilibrium 
from stochastic point of view. That is, the emergence of the phenotypic equilibrium 
is rooted in the stochastic nature of (almost) every sample path, the average-level 
stability just follows from it by averaging stochastic samples. 
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Figure 1: The phenotypic equilibrium of cancer cells. The figure is generated from the 
data (SW620 colon cancer cell line) in [|71. In this experiment, two cellular phenotypes 
were identified: cancer stem cells (CSCs) and non-stem cancer cells (NSCCs). It is shown 
that no matter where the initial state is (as four different cases shown in the figure), the 
CSCs proportion will converge to a fixed proportion as time passes. The same is true for 
NSCCs proportion. This phenomenon is termed phenotypic equilibrium [|6l. 


1 Introduction 

Stability is ubiquitous in biology, ranging from physicochemical homeostasis in cellular 
microenvironments to ecological constancy and resilience [[Il|2l|3. It is noteworthy that 
not only can the stability phenomenon arise in normal living systems, but it can also hap¬ 
pen in abnormal organisms such as cancer. As a large family of diseases with abnormal 
cell growth, cancer is generally acknowledged to be the malignant progression along with 
a series of stability-breaking changes {e.g. genomic instability) within the normal organ¬ 
isms a. However, some recent researches reveal the other side of cancer. An interesting 
phenotypic equilibrium was reported in some cancers [|5] |6l |71. That is, the population 
composed of different cancer cells will tend to a fixed equilibrium of phenotypic propor¬ 
tions over time regardless of initial states (Fig. 1). These findings provided new insights 
to the research of cancer heterogeneity. 


2 










The experimental works also stimulate theoretieians to put forward reasonable models 
for interpreting the phenotypie equilibrium (bl [8l IH [TOl HU [lU [T3l [Ml . In particular, it 
was reported that the intrinsic interconversion between different cellular phenotypes, also 
called phenotypic plasticity [fTSiri^ . could play a crucial role in stabilizing the mixture of 
phenotypic proportions in cancer. As a pioneering work, Gupta et al proposed a discrete¬ 
time Markov chain model to describe the phenotypic transitions in breast cancer cell lines 
dH. In their model, three phenotypes were identified: stem-like cells (S), basal cells (B) 
and luminal cells (L). The phenotypic transitions among them can be captured by the 
transition probability matrix as follows: 


/ 1 — Ps^B — Ps^L Ps^B Ps^L \ 

P = j Pb^s 1 — Pb^s — Pb^l Pb^l , (1) 

\ Pl^s Pl^b 1 — Pl^s ~ Pl^b ) 

where Pi^j represents the probability of the transition from phenotype i to j. Accord¬ 
ing to the limiting theory of discrete-time finite-state Markov chain, there exists unique 
equilibrium distribution vr = (tts, tt^, tt^) such that vr = vfP, provided P is irreducible 
and aperiodic [fTTl . The Markov chain will converge to if regardless of where it begins. 
By fitting the Markov chain model to their experimental data, the equilibrium propor¬ 
tions of stem-like, basal and luminal cells were predicted by the equilibrium distribution 
TiSi '^B, respectively. 

Even though the Markov chain model fitted the experimental results in breast cancer 
cell lines very well, Zapperi and La Porta dHll questioned the validity of the phenotypic 
transitions and gave an alternative explanation to the phenotypic equilibrium, which was 
based on the conventional cancer stem cell (CSC) model with imperfect CSC biomarkers. 
Moreover, Liu et al showed that the negative feedback mechanisms of non-linear growth 
kinetics of cancer cells can also control the balance between different cellular phenotypes 
[flSll . These works suggested that the phenotypic plasticity may not be the only explana¬ 
tion to the phenotypic equilibrium. To further reveal the mechanisms giving rise to the 
phenotypic equilibrium, it is more convincing to study the models integrating the phe¬ 
notypic plasticity with the other conventional cellular processes of cancer. Motivated by 
this, a series of works discussed the phenotypic equilibrium by establishing the models 
coordinating with both hierarchical cancer stem cell paradigm and phenotypic plasticity 
iiiioiiniiniiiiiii. In these works, the phenotypic equilibria were intimately related 
to the stable steady-state behavior of the corresponding ordinary differential equations 
(ODEs) models. In other words, if one can model the dynamics of the phenotypic propor¬ 
tions as the following system of ODEs 


dx 

dt 


F{x), 


the unique stable fixed point x* (if exists) corresponds to the equilibrium proportions. 
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The aforementioned works have showed that the phenotypie equilibrium ean be ex¬ 
plained by different eoneepts of stabilities in different models. Thus a natural question 
is whether there exists a unified framework to harmonize the equilibrium distribution of 
the Markov chain model and the stable steady-state behavior of the ODEs model. In this 
study, we try to address this issue by establishing a multi-phenotype branching (MPB) 
model [fT9l . On one hand, our model integrates the phenotypic plasticity with the cellular 
processes (such as cell divisions) that have extensively been studied in cancer biology. On 
the other hand, the model is stochastic and closer to the reality with finite population size 
[I 20 II 2 TII . Based on this model, it is shown that the ODEs model can be derived by taking 
the expectation of our model. More specifically, the ODEs model is just the proportion 
equation of the MPB model. Besides, the Markov chain model is also shown to be closely 
related to our model. That is, the Kolmogorov forward equation of the continuous-time 
Markov chain model is a special case of the proportion equation provided that the divi¬ 
sion rates of stem-like, basal and luminal cells are the same. Interestingly, “same doubling 
time” of the three phenotypes was just observed in Gupta et a/’s experiment when they 
used the Markov chain model to explain the phenotypic equilibrium which is in line 
with our theoretical prediction. Moreover, our result also shows that one should be more 
cautious about the application of the Markov chain in modeling cell-state dynamics in 
larger time scales, since the Markov chain model takes no account of different capabili¬ 
ties of divisions by cancer stem cells and non-stem cancer cells. 

More importantly, by showing almost sure convergence of the MPB model, the sta- 
tionarity of the Markov chain model and the stability of the ODEs model can be unified as 
the average-level stability of our model. Note that the almost sure convergence indicates 
the path-wise stability of stochastic samples, providing a more profound explanation to 
the phenotypic equilibrium. In other words, the phenotypic equilibrium is actually rooted 
in the stochastic nature of (almost) every path sample; the average-level stability just fol¬ 
lows from it by averaging all the stochastic samples. Eurthermore, it is also shown that, 
not only can the model with phenotypic plasticity give rise to the path-wise convergence, 
but the conventional cancer stem cell model without phenotypic plasticity can also lead 
to the convergence under certain conditions. This echoes the works [l8l[l8l that the phe¬ 
notypic plasticity is not the only explanation to the phenotypic equilibrium. 

The paper is organized as follows. The model is presented in Section 2. Main results 
are shown in Section 3. Conclusions are in Section 4. 
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2 Model 


2.1 Assumptions 

In this section we give the assumptions of our model. Consider a population composed 
of different cancer cell phenotypes. For pure theoretical investigations, the number of the 
phenotypes can be any n in general [fT3ll^ . However, to better illustrate our theoretical 
results on the basis of more concrete biological background, enlightened by we here 
focus on the specific model consisting of three phenotypes: stem-like cells (S), basal cells 
(B) and luminal cells (L). The main assumptions are listed as follow: 

1. Stem-like cells can perform three types of divisions: symmetric division, symmet¬ 
ric differentiation and asymmetric division [|^ [24l l25l . That is, a stem-like cell can 
divide into two identical stem-like cells (symmetric division) or two identical differenti¬ 
ated cancer cells (symmetric differentiation; it can also divide into a stem-like cell and a 
differentiated cancer cell (asymmetric division). 

• symmetric division: S S-i-S; 

• symmetric differentiation: S B-i-B or S L-i-L; 

• asymmetric division: S S-i-B or S S-i-L. 

as is the division rate (or termed synthesis rate lITSll f. with the meaning that a stem-like 
cell will wait an exponential time with expectation as and then perform one particular 
type of division with probability P* (note that ^ Suppose the waiting time 

and the division strategy are independent to each other, then the product of 0:5 and Pi 
governs the reaction rate of the corresponding division type. 

2. For non-stem cancer cells, i.e. basal or luminal cells, we assume that not only can they 
undergo symmetric divisions with limited times, but they can also perform phenotypic 
conversions. To illustrate this, let us take B phenotype as an example. Suppose a newly- 
born B cell can divide at most m times. If we denote Pj as the B cell that has already 
divided i times, then we have the following hierarchical structure: 

• Bo^Bi-rBi; 


R R _i_R • 


B™, ^ 0. 
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as is the division rate, and is the death rate of B^- Moreover, assume that a B cell 
can convert into an S cell (termed de-differentiation [l26l l by phenotypic plasticity. Let 
the dedifferentiation rate of Bj be then we have 

• Bo^S; 


• B^^S. 

For simplicity, it is often assumed that Pbo = ••• = f^Bm llI3’ denoted as /Sb for short. 
Meanwhile, note that a B cell can also convert into an L cell [[ 6 l|. Since the biological 
mechanisms of the phenotypic conversions between different non-stem cancer cells are 
still poorly understood, for simplicity it is assumed that the phenotypic transitions be¬ 
tween B and L can only happen in the same hierarchical level. That is, supposing that a 
newly-born L cell can also divide at most m times, Lj is the L cell that has already divided 
i times, then we have 

• U; 

7 b is the transition rate. In fact, this assumption implies B —> L with constant rate jb 
overall, which is in line with the assumption in [j 6 |. For luminal cells, similarly, their 
cellular processes are shown as follows: 

• Lj Lj+i-i-Li+i (0 < i < m - 1); 

• ^ 0 . 

• Lj S (0 < i < m); 

• Lj Bj (0 < z < m). 

2.2 Multi-phenotypic branching (MPB) model 

Based on the cellular processes listed in the last subsection, we can model this cellu¬ 
lar system as a continuous-time Markov process in the discrete state space of cell num¬ 
bers (Chapter 11 in IITTII L If we let Xi be the cell number of S phenotype, X 2 = 
{X 2 ^\ X 2 ^\ be the vector representing the cell numbers of Bj cells, and X 3 = 

(Xg^^ ..., Xg™^)^ bc thc vcctor representing the cell numbers of Lj cells, then the 

dynamics of X = (Xi, X 2 , Xs)^ can be modeled as a multi-phenotype branching process 
flU. If we define Pr(T'; t) be the probability of X = a; at time t, according to the theory of 
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Chemical Master Equation (CME), the rate of ehange of Pr(x; t) is equal to the transitions 
into X minus the transitions out of it, i.e. 

^ ^ ]t) T^_^f/Pr(f; t), ( 2 ) 

x'i^x x'i^x 

where is the transition rate from x' to x and is the transition rate from x to x' 
(seej^for more details). 

In next seetion we will show that the ODEs model and the Markov ehain model ean be 
derived from our model. Eor convenience we term our multi-phenotype branching model 
the MPB model. 


3 Results 

3.1 Deterministic equations derived from the MPB model 

To relate our MPB model to the ODEs model, we consider the mean dynamics of the 
MPB model by averaging all the stochastic samples of it. 

Eet ) be the expectation of X, that is, for each component we define 

{Xi) ;= y^a:^Pr(T; t). 

X 

We multiply Xi on the both sides of Eq. 0. and then calculate the summation over all x 

Y = X] I X] Ta?/^^Pr(f' ;t) -Y ^x-^x'Pr(^; t) 

X X \x'^x x'^x 



Eor S cells: 

j( Y \ m 

= as (Pi - P 2 - Ps) (Xi) +(3bY +f^LY • ( 3 ) 

2=0 2=0 

Eor B cells: 

[ T = (2P2 + P4) (Xi) - {as + ^B + 1b) {xY + 7LXf; 

J d{x^2 ) ^ 2aB{xY^^) - (qb + /3b + 1b) (X^'^) + IlX^^'^ (4) 

i = 2aB{x)rY - (aB^ +Pb + 1b) (X^) + IlX^■ 
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For L cells: 

[ T = ( 2 P 3 + n) (Xi) - («z. + /3 l + 7l) (Xf) + ); 

< = 2aL{xt"'’) - (ftL + /9l + 1l) (X«) + 7 sX« (1 < z < m - 1 ); (5) 

i = 2«z.(xf +p, + 7l) (Xf^) + 7 bX^\ 


Then it is not diffieult to see that the dynamies of (X) can be eaptured by a system of 
linear ODEs, 


d{X) 

dt 


= G(X), 


where 


^ ~ [5'7]{2m+3)x(2m+3) — 


<^s{Pl — P2 — P'i) 





( 6 ) 


Pb Pl ••• 

Pl\ 


(oiB+PB+'yB) 0 

... ... 

0 

• (7) 

2a B —{aB+PB+lB) 

0 . 


. Q describes the eell 

number 

dynamics 


of eaeh phenotype at each hierarchical level. If we denote X 2 = YlT=o ^ 2 ^ "^3 ^ 

^ 3 '’ numbers of B and L phenotypes respectively, then it is often 

the dynamics of X* = (Xi, X 2 , Xs)^ that interests people. That is, 


'^ = as{Pi-P2- P3) (Xi) + I 3 b{X 2 ) + / 3 l(X 3 ); 

= as {2P2 + P4) (Xi) + (as - Pb - 1b){X2) + 1l{X^) - {ub + a_B„)); 
^'^ = as (2P3 + P3) (Xi) + 7 s(X 2 ) + (as - Pl - 1 l){X^) - {aL + aUiX^^^). 


( 8 ) 

M\ 


We can see that Eq. (js) is not linear of (X*), whieh also depends on (X^™^) and (X 3 ™'’) 
separately. Teehnically this is due to the limited capability of divisions of B and E pheno¬ 
types. In the limit of m, or when m is relatively large in eomparison to observational time 
scales (e.g. t ^ m), Eq. M can approximately be expressed as a linear system of (X*): 


d{X*) 

dt 


where 




G*(X*), 


(9) 

/ as{P\—P2—P3) Pb 

Pl \ 


3 x3 = ( as{2P2+P4) aB-PB-lB 

IL I . 

( 10 ) 

V as{2P2,+P5) IB 

aL-Ph-lL J 



In this way the model reduees to the three-phenotypic model investigated in ifTTll . How¬ 
ever, Eq. ([^ should be adopted for describing larger time scales (e.g. t 3> m). Note that 
it is inconvenient to analyze Eq. ([^ direetly, we will show later that analyzing Eq. is 
quite helpful for the understanding of Eq. ([^, especially in the study of the phenotypic 
equilibrium. 
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3.2 Proportion equation: Bridging the MPB model and the ODEs 
model 

Since Eq. ([^ describes the dynamics of the absolute numbers of different eellular phe¬ 
notypes, we term it the number equation. However, to investigate the phenotypie equilib¬ 
rium, we are more concerned about the dynamies of the relative numbers {i.e. proportions) 
of different cellular phenotypes. Let p be the vector representing the proportions of dif¬ 
ferent cellular phenotypes. By replacing (X) in Eq. ([^ with p , we have the equation 
governing the phenotypie proportions as follows (see|^ 

^=Gp-pe^Gp, ( 11 ) 

where e = (1,...,!)^. We term Eq. (|TT ) the proportion equation. It is noteworthy that the 
stable steady-state behavior of Eq. ( [TT] ) just eorresponds to the phenotypie equilibrium 
investigated in ifT^ [T3l . The proportion equation thus eonnects the MPB model and the 
ODEs model in previous literature, implying that the ODEs model can be seen as the 
average-level eounterpart of the stoehastie MPB model. To show the stability of Eq. d). 
we have the following theorem (seefor the proof): 

Theorem 1. There exists unique positive stable fixed point jl in Eq. ( [77] ) provided that G 
is irreducible^ 

Theorem shows that the deterministie population dynamics of cancer eells will tend 
to an equilibrium mixture of phenotypie proportions as time passes. Besides, let p* be the 
proportion veetor of X*, i.e. 


f = iP*l,P*2,Pl) = 

Given lim4_).oo P = fi (Theorem[T]), 

mm mm 

= ^lnn(pi, = (/ii, 

2=0 7=0 7=0 7=0 

Thus we have the following result for fT: 

Corollary 1. Under the same condition in Theorem jT will tend to a fixed positive 
vector fi* as t ^ oo. 

' Strictly speaking, for completing the theorem it is necessary to add a small perturbation to the initial 
state in rare cases, see[C] 
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Corollary [^indicates the phenotypic equilibrium of the three-phenotypic model in Eq. 
Moreover, it should be pointed out that, the results in Theorem and Corollary can 
be seen as the average-level stabilities following from the the path-wise convergence of 
the MPB model, which will be discussed in Sec. 


3.3 The Markov chain model as a special case of the proportion equa¬ 
tion 

Note that the Markov chain model Eq. Q is discrete-time and the MPB model is continuous¬ 
time; to compare the two models in the same time scale, we turn our attention from 
discrete-time Markov chain to continuous-time Markov chain. Consider the standard 
model of continuous-time Markov chain. That is, let Pi{t) be the probability of the 
Markov chain being in state i at time t, its dynamics can be captured by the Kolmogorov 
forward equation: 


dP{t) 

dt 


Q^m, 


where Q-matrix [gjjjsxs satisfying 


( 12 ) 


Qij > 0 Vz 7^ j, 


(13) 


Qii ^ ^ Qij ■ 


(14) 


P 


We now discuss the relation between P{t) and p*. By replacing {X*) in Eq. (|9|) with 
', we obtain the proportion equation governing p* 


— = G p -p e G p , 


(15) 


where e = (1,1,1)^ and G* in Eq. (10). If we let the sum of each column of G* is the 
same, i.e. 


as = aB = as = zc, 


then Eq. ( [T5] ) becomes 

^The derivation of Eq. (15 1 is similar to that of Eq. (Ill, see[B| 


( 16 ) 
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where / is identity matrix. If we denote H = (G* — it ean be shown that H satisfies 
the eonditions (13) and ( [T4] ) for the Q-matrix (see|^. In other words, the Kolmogorov 
forward equation Eq. ([T^ is a special linear case of the nonlinear proportion equation 
Eq. l(igi. This relation implies that, when the division rates of the three phenotypes are 
the same, the dynamics of the phenotypic proportion can equivalently be captured by the 
Markov chain model where only the phenotypic transitions are accounted for. Otherwise, 
the Markov chain model may oversimplify the phenotypic dynamics with unequal division 
rates. Interestingly, it was reported in Gupta et a/’s experiment that the subpopulations of 
S, B and L phenotypes have the same “doubling time” dH, which justified their application 
of the Markov chain model. However, as mentioned in the end of Sec. 3.T Eq. S 
is valid only in relatively short time scale. Eor larger time scales, it is unreasonable to 
model the three-phenotypic dynamics by the Markov chain model taking no account of 
different capabilities of divisions by cancer stem cells (unlimited) and non-stem cancer 
cells (limited), even if they have the same division rate. Therefore, one should be cautious 
about the application of the Markov chain in modeling cell-state dynamics. 


3.4 Path-wise convergence of the MPB model 

We have seen that the MPB model provides a unified framework for the ODEs model and 
Markov chain model. In this section, we will show path-wise convergence of the MPB 
model, which provides a much stronger concept of stability by which both the stable 
steady-state behavior of the ODEs model and the equilibrium distribution of the Markov 
chain model will serve as average-level stabilities of the MPB model. 

Much attention has long been paid to the limit theorems of multi-type branching pro¬ 
cesses by mathematicians [l2^l29l[30l[3T]| . Here we are not going to discuss the rigorous 
mathematical theory in general (which is the focus of our another work [l22l k Instead we 
are more interested in the specific results related to the phenotypic equilibrium, i.e. the 
conditions under which p converges to a positive vector jl. Unlike the p in Theorem[^ the 
p here is stochastic. The “convergence” here means almost sure convergence. That is, if 
the convergence of p holds, almost all the stochastic paths will tend to a fixed equilibrium 
(also termed path-wise convergence). 

We present our main results in the following two theorems (see for the proofs and 
mathematical details): 

Theorem 2 . If G in Eq. ([^ is irreducible and its Perron-Frobenius eigenvalue is positive, 
then p will tend to a fixed positive vector fi almost surely as t ^ oo conditioned on 
non-extinction of the population. 

Theorem 3. Assume that 

(1) all the phenotypic transition rates are zero, i.e. fis, fii eif^d Jb, 'Jl eire zero; 
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( 2 ) as > 0 , as > 0 and as > 0; 

(3) Pi > 0 (1 < i < 5) and Pi > P 2 + P 3 ; 

then p wUl tend to a fixed positive vector fi almost surely as t ^ 00 conditioned on 
non-extinction of stem like cells. 

The above two theorems are applieable to different cases. Theorem corresponds 
to the case with phenotypic plasticity, since the irreducibility of G is satisfied as long as 
the conversions between different phenotypes can happen. In contrast, Theorem [^corre¬ 
sponds to the case without phenotypic plasticity, since all the phenotypic transition rates 
are assumed to be zero. Interestingly, even though the assumptions of the two theorems 
are basically different, both of them can lead to the path-wise convergence p. Further¬ 
more, it is easy to see that the path-wise convergence of p* is implied by Theorems [^ and 

s 

Corollary 2. p* will tend to a fixed positive vector fi* almost surely as t ^ 00 under the 
conditions in either Theorem^or Theorem\^ 

Figs. 2 and 3 illustrate the path-wise convergence of p* implied by Theorems [^ and [^ 
respectively by using stochastic simulations (j^ shows the simulations for p in details). In 
both cases, even though all the stochastic paths fluctuate at the beginning of the process, 
the proportions of S, B and L cells eventually converge to their equilibrium proportions 
as time passes. Since the path-wise convergence indicates the stability of (almost) every 
stochastic sample, the convergence of the mean dynamics just follows from it by aver¬ 
aging all the stochastic samples (see lower panels of Figs. 2 and 3). Note that both the 
Kolmogorov forward equation of the Markov chain and the ODEs model can be seen as 
the mean dynamics of the phenotypic proportions; their stabilities just correspond to the 
average-level stabilities of the MPB model, which can be seen as direct results of the 
path-wise convergence. In this way, the path-wise convergence provides a deeper under¬ 
standing to the phenotypic equilibrium from the stochastic point of view. 

As the end of this section, it is noteworthy to emphasize that, according to Theoremj^ 
the phenotypic equilibrium can still happen in the paradigm of conventional cancer stem 
cell theory. The assumptions in Theorem together indicate the cellular hierarchy pro¬ 
posed by the cancer stem cell theory [|24ll . That is, cancer stem cells (S cells) are capable 
of differentiation into other more committed non-stem cancer cells (B and L cells) but 
not vice versa. In this way, cancer stem cells are at the apex of this cellular hierarchy. 
Moreover, the assumption “Pi > P 2 + P3” implies the dominance of S phenotype during 
the growth of the population. To show this, note that (Pi — P2 — P3) is the eigenvalue 
corresponding to S phenotype of G in Eq. ( [T9| ), which is the only positive eigenvalue of 
G provided “Pi > P2 -f P3” (seej^. In other words, instead of the phenotypic plastic¬ 
ity, Theorem [^ also gives an alternative explanation to the phenotypic equilibrium in the 
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Proportion of each cellular phenotype in path 



Proportion of each cellular phenotype in average 



Time 


Figure 2: Stochastic simulations for the case with phenotypic plasticity (Theorem [^. 
Upper panel shows the stochastic path-wise dynamics of the phenotypic proportions of S 
(blue), B (black) and L (red). The initial numbers of S, B and L cells are assumed to be 20, 
0 and 0 respectively, that is, the initial proportions of S, B and L cells are 100%, 0% and 
0%. According to the assumptions in Theorem]^ we set m = 10; as = 0.8, Pi = 0.3, 
P 2 = 0.2, P 3 = 0.2, P 4 = 0.15, P 5 = 0.15; as = 0.6, as^ = 0.3, /3b = 0.1, 75 = 0.05; 
ai = 0.7, = 0.3, I3l = 0.13, 7 l = 0.2. Thirty stochastic samples for each phenotype 

were produced. It is shown that even though the stochastic paths fluctuate at the beginning 
of the process, the proportions of S, B and L phenotypes eventually path-wisely tend to 
their equilibrium proportions respectively. Lower panel shows the mean dynamics of the 
phenotypic proportions by averaging all the thirty samples shown in upper panel. 
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Proportion of each cellular phenotype in path 




Figure 3: Stochastic simulations for the case without phenotypic plasticity (Theorem [^. 
The initial cell numbers of S, B and L cells are also assumed to be 20, 0 and 0 respectively. 
According to the assumptions in Theorem we set m = 10; as = 0.8, Pi = 0.5, 
P 2 = 0.14, P 3 = 0.16, P 4 = 0.1, P 5 = 0.1; as = 0.4, = 0.3, 13b = 0, 75 = 0; 

ai = 0.45, ai^ = 0.3, I3l = 0, 7 l = 0. Ten stochastic samples for each phenotype were 
produced. Upper panel shows the path-wise convergence of the phenotypic proportions. 
Lower panel shows the average-level stability of the mean dynamics. 
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framework of the cancer stem cell theory, as long as the cancer stem cell phenotype is 
dominant in the population. However, it is interesting to see that the convergence rate of 
the case in Fig. 2 is faster than that of the case in Fig. 3, even though they both give 
rise to the path-wise convergence. This suggests that perhaps the convergence rate (rather 
than the convergence itself) could serve as an indicator to distinguish the models with 
and without phenotypic plasticity, which might be another meaningful research topic in 
future. 


4 Conclusions 

In this study, we have presented a multi-phenotype branching model of cancer cells. On 
one hand, this model can serve as an underlying model from which the ODEs model 
and the Markov chain model can be deduced. On the other hand, the almost sure con¬ 
vergence of the model enhances our understanding of the phenotypic equilibrium, from 
average-level stability to path-wise convergence. Furthermore, our results have indicated 
that, even though the phenotypic plasticity facilitates the phenotypic equilibrium, it is not 
indispensable in some cases. It has been shown that the conventional cancer stem cell 
model can also stabilize the mixture of the phenotypic proportions, providing an alterna¬ 
tive explanation to the phenotypic equilibrium. 

Moreover, it should be noted that even though this work is focused on the issue of 
cancer, our methods can conveniently be used to more generalized cell population dy¬ 
namics ^T2\ . To further reveal the biological mechanisms of the phenotypic equilibrium, 
more detailed dynamic models of cancer cells are needed. For instance, the hypothesis of 
cooperation among cancer cells has been put forward [[32l . In particular, self-sufficiency 
of certain growth signals of cancer cells supports the concept of mutualism and could 
be an important mechanism supporting the phenotypic equilibrium. Therefore, the mod¬ 
els of capturing the interactions among cancer cells, e.g. evolutionary game models (331, 
could be a promising research direction in future. Furthermore, the genetic and epigenetic 
state networks (341 l35l of cancer will enable us to explore the molecular mechanisms of 
the phenotypic equilibrium, which are poorly understood. The network methods have 
successfully been used to investigate the processes of cellular pluripotent reprogramming 
(3^ and epithelial-mesenchymal transitions (EMT) (T7l . Note that EMT could play a key 
role in regulating the phenotypic heterogeneity in cancer (381 . further studies on it should 
be another important tasks in future plans. 
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A Expanded form of Eq. Q 

Here we show more details about the master equation Eq. Q 

(iPr(T; t) 


dt 


T5j/^^Pr(f'; t)-^ t). 

x'^x 


To obtain the expanded form of Eq. Q, we need to confirm all possible and 

T^'^x- Based on the model assumptions, we can calculate and Tg/^g correspond¬ 


ingly. For example, let x' = 

±,^2 ,..., X 2 ,...,^3 ) 


as Pi 


[Xi - 1,X^2 
V r ( 0 ) 

-)■ (xi,a;2 S 


,..., X 2 
(m) 

...,^2 , 


,xf\...,x^^^)'^, the event “(xi — 
X 3 °\ ...,X 3 ™^^)^” will happen if any 


one of “S' —> S -f S'” happens. Since the number of S cells is xi — 1 in current pop¬ 
ulation, the transition rate should be (xi — 1) x asPi- On the other hand, the reaction 




(xi + l,x^°\ 


S + S'” can also lead to the transition from (xi, X 2 


(0) 


, X, 


M ^(0) 


.., x^^\xf\ ..., Xg”^'*)^ with rate Xi x asPi- Along this way we can deter- 


d^)\T 




.,X3 


(m)\T 


-)■ 
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mine all the transition rates similarly. Therefore, 

t) = 

x'^x 

{xi - l)asPiPr{xi - ...,4™^) 

+ (a;i + l)asP2Pr{xi + l,a;® - 2, ...,x^^\xf\ ...,x^^^]t) 

+ (a;i + l)asP3Pr(a;i + 1,•••, 4°^ - 2, ...,4'"4) 

+xiasPi^x{xi, x^^'> - 1,..., x^^\xf\ ...,4'"^; t) 

+xiasP:i^x{xi, 4°\ •••, 4™\ 4°^ - •••,4'"^; 

m 

+ + l)aBPr(a:i,...,4*”^^ + 1) ^ 2 ^ ~ 2,...;t) 

i=l 

+ (4™^ + l)ai?™Pr(a^i,4°\--,4™^ + 1,4°^ •••> 4™4) 

m 

+ y^(4*^ + l)/^BPr(a^i - l,---,4*^ + 


i=0 


+ X](4*^ + l)7BPr(a^i, •••,4*^ + •••> 4*^ - 1, •••; ^) 

m 

+ y](4*~^^ + l)ttLPr(a^i,--,4*~^^ + 1>4*^ - 2,...;t) 
+ ( 4 ”"^ + l)«L™Pr(a^i, 4 °\ •••, 4 ”"^ 4 °^ •••, 4 '"^ + 1 ; 

m 

+ y](4*^ + l)/5LPr(a:i - 1,..., 4*^ + 1,...; t) 


i=0 


+ y^( 4 *^ + l)7LPr(a^i,--,4*^ - l,--,4*^ + 


(7 


77 


i=0 


and 


5 

Tg^g^Pr{x-, t) =xi^ asPiPr{x-, t) 

x'^X 2=1 

m—1 m 

+ y] aBa:2*^Pr(f; t) + aij„a;4^Pr(T; t) + y (/^s + 7B)4^Pr(p t) 

i=0 2=0 


m—1 m 

+ y aLX^3^Pv{x-, t) + aL„a;4^Pr(f; t) + y (/^l + 7L)a^?Pr(p t). 
2=0 2=0 
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B Proportion equation and Kolmogorov forward equa¬ 
tion 


Firstly we show how to derive the proportion equation ( [TT] ) from the number equation Q. 
Let us rewrite the matrix form of Eq. (|^ into eomponent form|^ 


d{X,) 

dt 


gil{Xi) + gil{X2) + ... + gin{Xn). 


Note that 

{X.) {X.) 

(Xi+X2 + ...+X„) N ’ 

d{X,) _ d{p,N) _ dN d^ 

dt dt dt dt ’ 

then 


dpi 1 d{Xi) Pi dN 
dt N dt N dt 

n n n 

= J2^gijPj) -p^Yl J2^9kjPj)- 

j=l k=l j=l 


When turning the above eomponent form baek to the matrix form, we get Eq. 0 

dp -,Tn~ 

dt 

Similarly, we ean also get the proportion equation Eq. ([^ for p*" 

(] 

-*f< -*|: T/O* /^n\ 

—— = G p — p e G p . (17) 

dt 


In what follows we show how the proportion equation Eq. ( [T7] ) relates to the Kolmogorov 
forward equation of eontinuous-time Markov ehain. If the eolumn sums of G* are the 
same and equal to k, i.e. 

Oig — — OiL — ki, 


then 


df 

dt 


G*f -Kf = H^f. 


^The dimension of G is 2m + 3, for simplicity we let 2m + 3 = n. 


18 










where H = (G* — kIY' and / is identity matrix. For any i ^ j, it is easy to see that 


hij — 


9 


Note that all the off-diagonal elements of G* are non-negative, so h^j > 0, 
satisfying the condition Eq. (131. Meanwhile, ha = g*^ — k. Note that k = J2i 9i 


ji' 




ha = g* 


K = — 




9 ji 




satisfying the condition Eq. Hence H corresponds to the Q-matrix of a continuous¬ 
time Markov chain. 


C Proof of Theorem [T| 


First of all, we have two remarks on G in Eq. Q: 

• Note that the off-diagonal elements of G are all non-negative, we call G an ME- 
matrix (see Chapter 2 in [fTTl l. For sufficiently large r, G + tI is a non-negative ma¬ 
trix (I is an identity matrix). In other words, ME-matrix is essentially non-negative. 
The ME-matrix G is said to be irreducible if G -\- tI is irreducible 0 


• When G is irreducible, from Theorem 2.6 in ifTTll . there exists a Perron-Frobenius 
eigenvalue Ai satisfying that i) Ai is real and Ai > ReA for any eigenvalue A 7 ^ Ai; 
2) Ai is simple, i.e. a simple root of the characteristic equation of G\ 3) Ai is 
associated with (up to constant multiples) unique positive right eigenvector fl. Here 
we assume that jl is the normalized right eigenvector of Ai, that is, /ii -f + ••• + 

l^n !• 

Since Ai is simple, the solution of Eq. Q can be expressed as 





Alt 






j =2 1=1 


i=l 


(18) 


where Ai, A 2 , • • ■ A^ are the different eigenvalues of G, nij is the algebraic multiplicity of 

Xj, rj- is the corresponding eigenvector of Xj, Cj^i is determined by initial states. Suppose 
ci^i 7 ^ 0, since ReA* < Ai (i 7 ^ 1), 


t^+oo 


m '^j 

= 11+ lim yy 

t^+00 ^ ^ Cl 1 ^ ’ 

7=2 1=1 i=i 




= /i. 


non-negative matrix M is said to be irreducible, if for every pair of indices i and j, there exists a 
natural number k such that is larger than 0. 
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Thus 


P = 


{X) 


W/ci,ie 


Alt 




/i 


{Xi-\-X2 + Xn) (Xi + ... + + ...+//. 


= /i > 0. 


Before eompleting the proof, we need to diseuss the ease ci,i = 0. In this ease, the 
above argument does not work. However, sinee fluctuations are inevitable in real world. 
Cl 1 = 0 will hardly happen in reality. To show this, let f = 0 in Eq. ([T8]) 


j=2 1=1 


This is a linear equation of CjBy Cramer’s Rule we have 

_ det|5*| 

“ det|5| ’ 

where B = [/W ^ ^ ^], B* is just B with its first column replaced by (Xq). It 

is easy to add a small perturbation ev to (Xq), so that all the columns of B* are linear 
independent, hence ci^ 7 ^ 0 holds. 


D Proofs of Theorems |2| and El 

The proofs of Theorems and are both on the basis of the following lemma 

Lemma 1 (Theorem 5 in [|22l). Assume that the Perron-Frobenius eigenvalue Ai of G 
in Eq. 0 is simple and positive. Conditioned on essential non-extinction, p will tend 
to fl almost surely as t ^ 00 . p is the normalized right eigenvector of Ai, which is 
non-negative. 

For proving Theorems and firstly we need to explain the concept of essential 
non-extinction. We are not going to discuss the general mathematical definition of it (see 
Sec. 4.2 in ll^ l. In our MPB model, essential non-extinction specifically means non¬ 
extinction of the phenotype corresponding to the Perron-Frobenius eigenvalue Ai. For 
Theorem the assumptions (2) and (3) implies that the Perron-Frobenius eigenvalue of 
G is Qii = as{Pi — P2 — P3) (we will show this later). In other words, the essential 
non-extinction here just means non-extinction of stem-like phenotype. For Theorem 

^As far as we know. Theorem 5 in ll22l requires minimal constraint to the path-wise convergence of our 
concern. However, it should be noted that technically our main results can also be proved based on Theorem 
3.1 in Ea. 
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since G is irreducible, it is possible for any two phenotypes to (directly or indirectly) 
inter-convert into each other. In this case, non-extinction of one particular phenotype is 
equivalent to non-extinction of any phenotype. This implies that, no matter which phe¬ 
notype corresponds to Ai, to guarantee essential non-extinction, it is sufficient to assume 
non-extinction of the population in general. Therefore, the conditions provided in Theo¬ 
rems [2] or [3] ensure the essential non-extinction of the model. 

We now start to prove the two theorems. On one hand, we need to show that Ai 
of G is simple and positive in both theorems. On the other hand, since Lemma only 
concludes the non-negativity of fl, we need to further show the positivity of fl. The proof 
for Theorem]^ is straightforward, since we assume that G is irreducible, according to the 
second remark in|^ Ai is simple and /i is positive. Note that Ai is also assumed positive, 
by Lemma [T] we have p ^ jl > 0 almost surely as f —)■ cx). 

For Theoremaccording to the assumptions, G reduces to a lower triangular matrix 
as follows 


G - [gij] 


/ as(Pi-P 2 -P'i) 0 . 

as{‘iP2+Pi) -oiB 0 ■■■ 

0 2q5 —olb 0 


as(2P3+P5) 0 

V .0 ::: ::: 


-OtL ■■■ 
2a.B —OIL 


0 

0 

0 

0 


\ 

/ 
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It is easy to know that the eigenvalues of G correspond to the diagonal elements. By 
assumptions (2) and (3), Ai = as {Pi — P2 — P3) is the Perron-Frobenious eigenvalue 
which is positive and simple. By Lemma we have p ^ fl almost surely as f —)■ cxd, 
where fl is the normalized right eigenvector of Ai. To complete the proof, we need to 
show that fl is positive. Note that fl satisfies the following equation 


G pL — Ai/i. 


( 20 ) 


By expanding this equation, we have 

Ai/ii = Ai/xi; 

(2F2 -f P4) /ii = {aB + Ai)/X2; 

2aBP2 = (tts + Ai)/r3; 

< 2aBPm = (21) 

as (2P3 -f P5) pi = {as + \i)p{m -\- 2 ); 

2ai,Pm+2 {c^L P Ai)/im-|-3) 

^2aLP2m+2 = + ^l)l^2m+3- 
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Suppose fii > 0, then we have 


/^2 


Ois (2-P2 + Pi) 
+ -^1 


fll>0 


sinee as ( 2 P 2 + Pi) > 0 and as + Xi > 0. With the same logie, we ean show the 
positivity of fii recursively, which completes the final proof. 


E Stochastic simulations for Theorems |2| and SI 
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Figure 4: Illustration of Theorem The parameters are the same as those in Fig. 2. 
Dynamies of B and L phenotypes at eaeh hierarehieal level are shows in (a) and (b) re- 
speetively. 
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Figure 5: Illustration of Theorem The parameters are the same as those in Fig. 3. 
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